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GENERATING UPWARD SWEEPS IN POPULATION USING THE 
TURCHIN-KOROTAYEV MODEL 

RICHARD E. NIEMEYER AND ROBERT G. NIEMEYER 


Abstract. The works of ICha-DunAlvInoNieCarFieLawllCha-Duril describe upward sweeps in populations 
of city-states and attempt to characterize such phenomenon. The model proposed in both ITurKorl ITurl 
describes how the population, state resources and internal conflict influence each other over time. We show 
that one can obtain an upward sweep in the population by altering particular parameters of the system of 
differential equations constituting the model given in ITurKorl ITurl . Moreover, we show that such a system 
has a nonstable critical point and propose an approach for determining bifurcation points in the parameter 
space for the model. 
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1. Introduction 

In |TurKor) . the authors attempt to construct a model for predicting the interrelationship between pop¬ 
ulation size (A(t)), state resources {S{t)) and internal conflict {W{t)). They build their model from first 
principles in Population Ecology and Economics, which are sound within those paradigms. We explain the 
necessary mathematics for an audience wishing to understand the Turchin-Korotayev models and incorporate 
similar techniques into their own research. More importantly, we examine the effects of changing various 
parameters in Equation © below. 
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where tq, po, a, a, 5, j3, c and 6 are parameters, whose interpretations are discussed in Sj3l As indicated 
above, the variables N, S and W stand for the population, state resources and internal conflict, respectively, 
and each is dependent on time t. The model given in Equation m is one that models the interaction between 
N, S, W of an agrarian society. 

The authors of [Turi ITurKor] do not list all of the parameters in their respective papers and do not explain 
how to force all of the variables to be nonzero, this being largely dependent on the technique one uses for 
numerically solving the system of equations. We give values for all of the parameters as told to us by P. 
Turchin via private communication. We provide Matlab code so that others may reproduce our results and 
those of [TurKor] . 

The models given in [Turl ITurKor) are meant to inspire those in the sociological sciences to find rigorous 
mathematical models for what are, to the trained eye, clearly dynamical systems. We attempt to convey the 
same necessity for formulating sociological problems as dynamical equations and modeling particular systems 
with more rigor as well as provide an explicit example providing a link between the work in [Cha-Dunl 
ICha-DunAlvInoNieCarFieLaw] on upward sweeps in historical data and the related work of P. Turchin and 
A. Koroteyev. 


2. Basic mathematical models of population growth and dynamics 

We want to refresh the reader on various concepts in this section on topics that are essential to the 
remainder of the article. We also want to demonstrate how it is one can model particular phenomenon in 
the sciences using differential equations. Both examples will become more relevant in SJ31 

Let us suppose we have box A and box B, each containing objects and there being no repeated objects 
in either box. A function is a rule for assigning an object in box A to one and only one object in box bQ 
For example, a rule for assigning a car to an owner would assign to each owner a car and such a car could 
be assigned to only one owner. Naturally, an owner could own more than one car, but only one person can 
be listed as the owner of a car. 

As an example, the mathematical rule that assigns to y one-third of the square of a number x is written 
as f{x) = and has the shape of what is called a parabola; see Figure [TJ Box A in this case is the set of 
real numbers (—oo, oo) and box B is the set of nonnegative real numbers [0, oo). 

To be technically correct, we refer to the mathematical phrase “y = /(a;)” as an equation. The derivative 
of a function, when it exists, is written as f'{x) and the equation y' = g{x) is called a differential equation. 
The differential equation y' = g{x) is called such for the fact that it is an equation involving a derivative. 
Now, sometimes, it is possible to determine the function f{x) such that f'[x) = g{x). But, for the most part, 
in practice it is quite difficult to calculate the antiderivative of a function, thereby solving the differential 
equation. 

This, however, does not prevent us from numerically solving many differential equations. There are a 
number of techniques for finding approximate solutions to a differential equation, many of them being built 
into the computational mathematics software package called Matlab. For some mathematical calculations, 
the software package R or Stata or SPSS are sufficient. If one wants to begin modeling certain social systems 
using advanced mathematics, it is necessary to use Matlab (or Octave, a freely available software package 
that, to put it succinctly, gets the job done and done well). 

An example of a differential equation that has an explicit solution is given by y' = 2x -I- 5. We can easily 
calculate that the antiderivative of 2a; -I- 5 is a;^ -I- 5a: -I- C, where C is an arbitrary constant of integration. If 
we suppose that yo = /(O) = 3, then (7 = 3. 

^There is such a notion of multi-valued functions, but we will never have a need to discuss those in this paper. 
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Figure 1. 


The graph of the function /(x) = \x^. We show only a portion of 


the graph of f{x). 


A more difficult to solve differential equation is y' = y. This is a function y = f{x) that is its own 
derivative. The only function that behaves in this way is the exponential function Ce^, where C is some 
constant. We can perform the following calculation. 


y 

(3) log y = 

glogy _ 

(4) y = 

This is made possible by the Fundamental Theorem of Calculus. A differential equation, like the one 
given at the beginning of Equation Q , has a general solution. A specific solution can only be given if we 
know how the function y(t) behaves at time t = to, the time from which the model is starting. The value of 
y at to is known as the initial condition. If j/(0) = 1 is the initial condition of the solution to Equation (P|). 
then y{t) = e*. 

It is reasonable to want to be able to specify an initial condition for a dynamical system. Otherwise, one 
is, in a sense, viewing all possible solutions at once. In Equation (|T]), all solutions to the differential equation 
given in Equation ([2]) are essentially listed simultaneously. The specific solution for when j/(0) = 1 is then 
the y{t) = e*, this being the solution to the initial value problem. 

Remark 1 (Variable vs. Parameter). A variable in a differential equation is a quantity changing with 
respect to time t or is the variable t itself. A parameter is some fixed quantity that is not changing with 
respect to time t. 

Example 1 (Variable vs. Parameter of a differential equation). In Equation ([5]), y is a variable. Since 
y = 1 ■ y, there is just one parameter in the differential equation, namely the coefficient of y. Implicitly, t is 
a variable since y is changing with respect to t. It may be possible for parameters to change, but not with 
respect to t. We will see this in the sequel. 

The graph of the solution to the initial value problem 


y 

1 

t + C' 

e^'e‘ 

Ce* 
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Figure 2. The graph of the function f{x) = e^. f{x) is the solution to the differential 
equation y'{t) = y{t). We have evaluated over the interval [—1,3]. 


Tn 

To 

Ti 

T2 

Ts 

T 4 

Ts 

To 

Tt 

Tg 

Tg 

Tio 

Tn 

Ti2 

Ti3 

Ti4 

Ti5 

Tie 

Ti7 

Tis 

Tig 

T 20 

T 21 

T 22 

T 23 

T24 

T 25 

R 

100 

95 

90 

81 

50 

32 

18 

11 

7 

13 

18 

19 

25 

40 

53 

71 

92 

105 

111 

125 

130 

138 

129 

115 

103 

98 

C 

100 

114 

120 

135 

143 

157 

147 

130 

120 

no 

107 

90 

79 

65 

59 

48 

34 

30 

36 

42 

53 

65 

83 

92 

99 

105 


Table 1. Recorded values for the populations of rabbits and coyotes in a particular environment. 


y'{t) = y{t) 

y(0) = 1 

is given by Figure [2| 

A Predator-Prey model is a standard example used to introduce the notion of a system of differential 
equations. It also serves as an example of our ability to capture key details about a seemingly complicated 
interaction between two species, namely a predator and its prey, using mathematics. Of course, the more 
equations you have involving the species, the more complicated the mathematics may become, but still 
mathematics provides explanatory power. 

Example 2. Consider two populations: rabbits and coyotes. Suppose an ecologist can record the increase 
and decrease in each population over a particular interval of time, say [1,7/] at equally spaced moments Ti, 
i = 1, ...,n, with Ti = 1 and r„ = Tf. The population data may appear as described in Table [H 

Upon examining the data, an ecologist may describe the interaction between rabbits and coyotes (which 
is clearly not a friendly interaction) in the following manner. 

We have been able to initially detect 100 rabbits (R) and 100 coyotes (C) in the region. They have 
all been tagged so that we may track their movement and determine at intervals of time the 
population of both rabbits and coyotes. We see from the data that there is a reciprocal relationship 
between the rabbits and the coyotes, the rabbits being prey for the coyotes. When the population of 
the coyotes increases to the point where the dwindling population of rabbits cannot sustain the larger 
population of coyotes, the coyotes begin to die off from lack of food and malnutrition. Consequently, 
the rabbit population eventually recovers and increases. With less coyotes to compete with, the 
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remaining population of coyotes can begin reproducing and picking off rabbits from the abundant 
population. 

Such a discursive description of the data is accurate, but can be summarized using mathematical equations. 
We can all agree that there is an interdependence between the two species. What we want to now demonstrate 
is how to model such a situation, in general, and how one may then search for parameters that recreate the 
situation illustrated by the data. We are not trying to fit an equation to the data. We are trying to condense 
the discursive explanation using mathematical equations. Computers understand equations, not discussions, 
and performing a mathematical analysis of the equations—that we will all agree upon shortly model the 
situation well—will provide insight into the dynamics of the system. 

In the absence of coyotes, for all intents and purposes, the rabbit population would increase without 
bound. Moreover, the rate of increase of the population can be reasonably modeled by assuming a rate of 
increase that is proportional to the population at the time of increaseQ The equation that relates the rate of 
change to the proportion of the population is given in Equation ([2|) and, hence, has a solution Ce*, for some 
constant C. Call this proportionality constant a. Of course, if the coyote population is not zero, then coyotes 
will hunt rabbits. For each rabbit, at time t there are C{t) many coyotes with which it could interact. Hence, 
there are RC many possible interactions. One can reasonably assume that with each possible interaction, a 
rabbit will either die or escape and live. Hence, some proportion of the number of interactions RC results 
in a rabbit’s death. This proportionality constant represents the efficiency with which a coyote can kill per 
unit of time. Denote this proportionality constant by (3. We can then say that the rate of change in the 
number of rabbits at time t is given by 


(5) — = aR- PRC. 

Now, the equation describing how the population of coyotes is changing at time t is not exactly similar 
since the coyotes are preying upon the rabbits. The rabbits need only reproduce and consume grasses and 
grains to ensure their survival (we are assuming an infinite supply of food). The coyotes’ food source is 
not infinite and restricted to a rabbit-only diet (assuming there are no other meaty animals around) and 
they must reproduce to increase their population. The rate at which the coyote population can reproduce 
is directly proportional to how much they can eat. Supposing the proportionality constant is 7 , the rate at 
which the coyote population increases is jCR. Now, the food supply for the coyotes is not endless and in the 
event there are no rabbits, we would expect the coyote population to decrease at a rate that is proportional 
to the population!^ Hence, we may write the proportionality constant describing the rate at which the coyote 
population is dying at time t as S and derive the equation below describing the rate at which the coyote 
population is changing. 


( 6 ) 

We may then combine Equations ([5]) and ([ 6 ]) into a system of differential equations that describes the 
coupled interaction of the two species. 


(7) 


dR 

dt 

dt 


= aR- PRC 


= -fCR - 6C. 


The system of differential equations given in Equation © is more generally known as the Lotka-Voltera 
equations and generally describes the interaction between a predator and a prey. In general, one can calculate 
a numerical solution to a system of nonlinear differential equations, modulo a few caveats. The Lotka-Voltera 
equations lend themselves well to being solve numerically. In Figure [3l we have plotted the population of 
rabbits against the population of coyotes at each time from Table[TJ One may then compare this to Figure 
S] illustrating the fluctuations of the populations over time under the explicit assumptions stated above. 


^This is called exponential growth and is seen in too many settings to count. 

^The coyote population is negatively impacted as a result of no food or fetuses die as a result of malnutrition. 
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Figure 3. This graph depicts how the population of the rabbits is related to the population 
of coyotes in our example. The data points are connected by line segments to illustrate the 
continuity in the population changes. 



Figure 4. The graph of the two solutions describing geometrically how the predator (coy¬ 
ote) and prey (rabbit) populations evolve in relation to the other. The x-axis is the prey 
population and the y-axis is the predator population. 


We want to examine now what happens when one makes adjustments to the parameters in the predatory- 
prey model in the example above. In particular, at t = 500, 7 changes from 0.1 to 0.2; at t = 1000, 7 changes 
from 0.2 to 0.3; at t = 3000, 7 changes back to 0.1. The phase portrait for this dynamical system is shown 
in Figure [HI 

In the graph describing the fluctuations in the rabbit population, we see that at t = 500,1000, 3000, the 
population of rabbits responds accordingly, as shown in Figure [ 6 ] We will argue in the following section 
that the upward sweeps in population described in |Cha-Dunl Cha-DunAlvInoNieCarFieLaw] are the result 
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Figure 5. When we change the parameters in the predator-prey model in a discrete fashion 
during the evolution of the system, we see that the populations follow different trajectories. 



(a) At t = 500, 7 changes from 0.1 (b) At t = 1000, 7 changes from (c) At t = 3000, 7 changes from 

to 0.2. 0.2 to 0.3. 0.3 to 0.1. 

Figure 6 . We adjusted 7 in the predator-prey model. These three graphs describe the 
fluctuations in the prey population subject to a change of the parameter 7 in the equation 
describing the rate at which the predator population is changing. The interval of time 
is [1,3500] Note, the scale of each graph is different so as to highlight the change in the 
solutions at the indicated time t. 


of a discrete change in particular parameters of the model described in the introduction and given in [Turl 
ITurKor) . 


3. The Turchin-Korotayev model, parameter changes and upward sweeps 

The system of differential equations given in Equation m with the initial conditions iV(l) = 1, S'(l) = 0 
and W{\) = 1 has a numerical solution shown in Figure |7| The parameter values are listed in Table HJ the 
duration of the model was 4,000 units of time. The interpretation of the parameters is given as follows. 

The parameter cq represents the population’s intrinsic growth rate, which is the maximum rate a popu¬ 
lation increases in size under ideal conditions (e.g. the population has access to unlimited resources, space, 
and perfectly ambient environmental conditions; there is an absence of predators, etc.). Broadly speaking, 
the ro is calculated in ecology as the population’s birth rate minus the population’s death rate. Accordingly, 
a population will grow if the birth rate is larger than the death rate; and a population will shrink if the 
death rate is larger than the birth rate. 
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0.015 
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0.05 

Po 

1.0 

/3 

0.25 

c 

2.0 

(5 

0.1 

a 

0.01 

a 

0.1 

^max 

3.0 




Table 2. The values of the parameters used in [TurKorl ITur 
eters prior to alteration. 


and the values of the param- 


Turchin and Korotayev [TurKor] theoretically derive tq from their integration of Thomas Malthus’s Prin¬ 
ciple of Human Population Growth and David Ricardo’s Theory of Marginal Returns; see [Mai] and [Ric] . 
respectively. According to Turchin and Korotayev’s reading of Malthus ( [Mai] 1. human populations grow 
exponentially as: 


( 8 ) 



wherein the per capita rate of population increase r is defined as a linear function of the per capita rate of 
surplus production, p{N), and a proportionality constant, C 2 : 


(9) r = C2p{N). 

Theoretically this means that the value of a population’s intrinsic growth rate is determined by the 
amount of surplus resources that each individual in the population is able to produce. In Turner’s ( [Trn ] ) 
reading of [Malj per capita surplus production is positively associated with technological innovation, such 
that an increase in available technology expands both a population’s productive capacity and level of available 
resources. Accordingly—least in this model—the intrinsic growth rate of a human population is significantly 
determined by technological innovation, such that an increase in innovation is positively associated with an 
increase in a population’s rate of growth. 

But, according to Turchin and Korotayev’s reading of Ricardo ( |Ric] b the per capita rate of surplus 
production p{N) is a declining function of population size, modeled as: 

( 10 ) p^N)= 

where k is the population size at which point the amount of surplus value produced equals zero, and ci 
is another proportionality constant. Theoretically this means that initial increases in population size will 
increase the rate of surplus production, but soon subsequent increases in population will lead to a decrease 
in the surplus production rate. Hence, Ricardo’s Theory of Marginal Returns qualifies Malthus’s Principle 
of Human Population Growth’s unrealistic assumption that the benefits of a single innovation continue 
indefinitely. 

Ultimately, in [TurKor] . Turchin and Korotayev integrate Ricardo’s and Malthus’s arguments by calcu¬ 
lating ro as the product of the two proportionality constants, ci and C 2 : 


(11) ro := C 1 C 2 

Thus, a positive increase in rg is theoretically interpreted as an increase in a population’s intrinsic growth rate 
due to an increase in the per capita production of surplus resources. This increase is enabled by technological 
innovation. 

The parameter po represents the state’s per capita taxation rate, calculated as the product of the pro¬ 
portionality constants ci and C3. As already defined above, ci is proportionality constant associated with 
Ricardo’s Theory of Marginal Returns. The proportionality constant C 3 is the proportion of surplus produc¬ 
tion collected by the state as taxes. Thus an increase in po reflects an increase in the state’s tax rate, and a 
decrease in po corresponds to a decrease in the tax rate. 
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1 < t < 1000 
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1000 < t < 2000 

5.0 

0.095 

0.45 

2000 < t < 4000 

7.0 

0.15 

0.95 


Table 3. The values of t for which /cmax, ro and <5 are changed. 


State resources, S - - Population, N Internal conflict, W| 



Figure 7. The solutions to the system of differential equations given by Equation (HD and 
parameters listed in Table [2] 


The parameter (3 represents the per capita state expenditure rate. Turchin and Korotayev assume the 
size of j3 is proportional to population size because an increase in population is generally associated with 
an increase in the amount of resources needing to be spent on the army and police, public works, and state 
bureaucracy. An increase (decrease) in /3 reflects an increase (decrease) in per capita state expenditures. 

The parameter a is a proportionality constant representing the frequency at which an encounter between 
to parties will result in violence. Likewise, the parameter b represents the rate at which members of each 
population are willing to ‘forgive and forget’ past grievances. The parameter b is proportional to war intensity 
because Turchin and Korotayev assume severe war leads to ‘violence fatigue,’ and thus a willingness to ‘bury 
the hatchet,’ so to speak. The parameter represents the effectiveness with which the state is able to suppress 
violence. 

Finally the parameters c and 8 represent the severity with which war affects the environment’s carrying 
capacity and the population size, respectively. The carrying capacity of the environment is given by fcmax in 
Equation CD- 


3.1. Upward sweeps in population from parameter changes. We wish to examine the effect of dis¬ 
cretely changing the parameters in the system of equations given by Equation (HD- We argue that one can 
recreate the phenomenon of the upward sweeps in population discussed in [Cha-DunllCha-DunAlvInoNieCarFieLaw] 
by changing particular parameter values in the system of equations given in Equation (HD. We examine the 
behavior of the solutions of an altered system of equations over the same time interval [1,4000] with the 
same initial values A^(l) = 1, S')!) = 0 and W{1) = 1. 

We see in Figure |8] that the population does exhibit an upward sweep in population when these three 
parameters are altered as described in Tabled 
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I state resources, S —Population, N — Internal conflict. W | 



Figure 8. The solutions to the system of differential equations given by Equation o and 
parameters altered as described in Table |31 


State resources, S - - Population, N Internal conflict, W| 



Figure 9. The solutions to the system of differential equations given by Equation ([T|). This 
time, we altered only Atmax as described in Tabled 


When we alter only the parameter fcmax, we see in Figure[9]that the population increases dramatically, buE 
relative to the graph of the numerical solution of W, there is not much of an upward sweep in populationjj 
Though, qualitatively speaking, one does see an upward sweep in the population. 

Altering only the the value tq, we do see in Figure [10] an upward sweep in the population at the times 
t = 1000 and t = 2000, as we have see injHJ but the intensity of internal conflict has not decreased. 

Finally, altering only S as described in Table jS] we see in Figure |Tl] that S and N are almost perfectly 
preserved, indicating that decreasing internal conflict does not have much of an affect on the population, 
thereby indicating no upward sweeps in the population, as one might expect. 


3.2. Alternate parameter changes for generating upward sweeps in N{t). We show that changing 
particular parameters in dS/dt and dW/dt at f = 1000,2000 results in an upward sweep in the population 
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We seek to characterize exactly what constitutes an upward sweep in population. 
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State resources, S — Population, N —Internal conflict, W| 



Figure 10. The solutions to the system of differential equations given by Equation ([T|). 
This time, we altered only rg as described in Tabled 


Slate resources, S — Population, N Internal conflict, W| 



Figure 11. The solutions to the system of differential equations given by Equation (IT|). 
This time, we altered only <5 as described in Table [S] 
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l<t< 1000 
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0.25 
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1000 < t < 2000 

5.0 

0.0833 

0.3 

2000 < t < 4000 

7.0 

0.0277 

O.T 


Table 4. The values of t for which fcmax, /? and po are changed. 


N{t) at the same times t. In particular, suppose we alter the parameters /cmax, P and po as described in 
Table m Then, Figure IT^ shows an upward sweep in the population at time t = 1000 and t = 2000. 

Additionally, changing the parameters of dW/dt as described in Table [5l we see that we again get an 
upward sweep in the population at the times at which the parameters are changed; see Figure [131 
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I state resources. S - - Population, N —Internal conflict, W| 



Figure 12. The solutions to the system of differential equations given by Equation o and 
parameters altered as described in Table IH 



a 
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a 

1 < t < 1000 

0.01 

0.05 

0.1 

1000 < t < 2000 

0.003 

0.15 

0.3 

2000 < t < 4000 

O.OOT 

0.45 

0.9 


Table 5. The values of t for which a, b and a are changed. 


I State resources. S -- Population, N —Internal conflict, W| 



Figure 13. The solutions to the system of differential equations given by Equation o and 
parameters altered as described in Table [SJ 


4. A STABILITY ANALYSIS 

4.1. The stability of Equation {[T]) with parameters given in Table For the parameter values given 
in Tabled! the system given in Equation o has the critical point 


( 12 ) 

(13) 


{N, S,T) 


p51 118401 3 \ 
Viw’ 256000’ 
(2.1938,0.4625,0.0375) 


The Jacobian of the system in Equation o with the same parameters as above, when evaluated at the 
critical point, is 
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(14) 


-0.2362 \ 
-1.1250 
-0.0500 / 

The eigenvalues of the Jacobian are approximately 


9 

0 

189 \ 


f - 0.0112 

_T 

0 

-T 

l-l 

-0.7500 

ssf 

8000 

1 

10 

-4 J 

20 / 

1 1 

0.0439 


Ai Rs -0.4085 

A 2 « 0.1736-0.10071 

A 3 Rs 0.1736+ 0.1007i 

In other words, the real part of the one of the eigenvalues is positive, meaning that the critical point is not 
a stable equilibrium. This means for any initial condition (A^(l), S')!), IT(1)) close to the equilibrium point, 
the trajectory of the solution will diverge from the equilibrium point. However, the solution approaches a 
limit cycle. In Figure 0 we see that the solutions are periodic and, were we to display the figure, one would 
see that the system is approaching a limit cycle in 3-space (with (TV, S, IF) at time t being plotted along the 
cc, y, z-axes). We refrain from showing the phase diagram of the solutions, because they are not as visually 
illuminating as the plots of the solutions themselves. 


4.2. Bifurcations in phase space. In H4.I1 we saw that Equation ([T]) with parameters given in Table [2] 
had an unstable equilibrium solution and that for any initial condition near the critical point, the graph of 
the solutions in phase space approached a limit cycle. 

Solving for the critical point {Nc, Sc, Wc) in terms of the parameters yields the following convoluted 
expression 


(15) 




c/3 

Po^ 





Po<5 


lVc = 


PoS 


that is valid under the following condition^ 


(16) 

(17) 


l35kcna.yiP0 + /3croPo < Cro/3^ -I- 

2al3^c^ro^po + 2a0^c5k^s.^ropQ + 2al3c5kcna.^r^P[)^ + 2a/3(5^fcmax^Po^+&/3<5ropo^ 
(Xj3 C Tq Clf3 C Tq Pq 4ft/3 cJ/Cmax^OPO 4” Clf3 5 TCmax Pq 


< 

-\- ClS kcciax Po . 


The characteristic equation for the Jacobian matrix is: 


are also assuming the parameters are positive. 
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4a/3^c^Aro^ — lOa^^c^Aro^po “ 6a/3^ci5fcniaxAropo “ Q:f3^cSroPo^ + 8a/3^c^Aro^po^ + 14a/3^ci5A:maxAroPo^ 

^ -aP'^cSXpo^ + 2aP'^cSropo^ + 2a^2(5^A:max^Apo^ + a(3‘^5'^k^a.^pQ^ - 2aPc^Xro"^po^ 

Spt 

-10a/3c5fcniaxAropo^ + 2a/3c5Apo^ - af3c6ropo'^ - 4a/3(5^fcniax^Apo^ - 2aPS'^kmaxPo'^ 

(-L®) + r 4 

^Po 

-PSX^roPo^ - bpSXroPo^ + 2ac5k^^^Xropo^ - acSXpo^ + 2a(5^A:max^Apo'‘ + a^^fcmaxPo® 

Spt 

SX^Po^ + SX'^ropo'^ + bSX'^po'^ + bSXropo'^ 

^ w 

As one can see, trying to compute the eigenvalues of the Jacobian matrix is no easy task. Instead, one 
should investigate whether or not a reduction in the number of parameters is possible. But, if one can 
determine that such a polynomial always has a positive root or a complex root with a positive real part, 
subject to the conditions stated in lines [16] and 1171 then we know that the system will always converge to a 
limit cycle and never have any bifurcation points in the parameter space. 

5. Discussion 

We have successfully demonstrated a sufficient condition for an upward sweep in the population by manip¬ 
ulating the intrinsic growth rate tq , carrying capacity for a system fcmax and the efficiency with which internal 
conflict W can affect the rate at which the population is changing. At first, increasing the efficiency with 
which internal conflict can affect the rate of population change may seem counter-productive for producing 
an upward sweep in the population. However, the change in the carrying capacity apparently offsets any 
negative consequences of this for the population and also results in the city-state governing more efficiently 
(a decrease in S(t)) and internal conflict being dramatically reduced, as well. Granted, the frequency with 
which population and internal conflict change is more dramatic, but the overall changes are less dramatic. 

The next step in this research program is to use the Turchin-Korotayev model for agrarian societies to 
simulate the rise and fall of two agrarian societies interconnected by trade T. Really, the variable T could 
represent any quantity shared or transferred by the two societies. We begin our next analysis by assuming 
T is trade and there are two nodes in the graph with an edge connecting them allowing for bidirectional 
transport of T. 


6 . Matlab code 

6.1. Predator-prey model with parameter adjustment. 

7o7oThese are the differential equations and paramieters. 
7o7oRobert Garrett Niemeyer 
7o7oRichard Evan Niemeyer 

function dy=predprey(t,y) 


7o7oThese are the parameters, as listed in your Figure 7 

alpha=5; 

beta=.1; 

gcmima=. 1; 

delta=5; 


if t>=50 && t<=100 
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gcumna = . 2; 

elseif t> 100 && t<=300 
gcunma = . 3; 
elseif t> 300 
gcunma = . 1; 


end 

dy=zeros(2,1); 

dy(l) = alpha*y(l)-beta*y(l)*y(2); 
dy(2) = gamma*y(2)*y(l) - delta*y(2); 


end 


7,7,solve_problem will run and produce three plots 

7.7.Syntax of solve_problem: solvePredPrey(lNlTlAL_PREY_POP, 1N1T1AL_PRED_P0P) 

7o7oRobert Garrett Niemeyer 

7o7oRichard Evan Niemeyer 

function [T,Y]=solvePredPrey(initCond) 

7o7oWe set the error tolerances very low 
options = odeset(’RelTol’,le-5,’AbsTol’,le-10); 


[T,Y] = ode45(@predprey,[1,350],initCond,options); 
7o7oWe plot the solutions 


plot(Y(;,1),Y(:,2)); 


end 

6.2. Equation {[1]) with parameter changes. The following three functions typed in Matlab code should 

be placed into their own files so-named. For example, the 

runParams 

function should be placed into a file called 
runParams.m 


7o7oRobert Garrett Niemeyer 
7o7oRichard Evan Niemeyer 
function runParams(initcond,params) 
format long 

7oWe want to know what parameters are going to be 
7»chariged at particular times and which are not 
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global paramsToTurnOnk; 

“/oWe generate all permutations of the paramters 
7,in the event one wants to examine 343 different 
7oWays to alter the parameters at particular times 

paramArray = [1,2,3,4,5,6,7,8,9]; 
paraunsToTurnOn = zeros(343,9); 
paramsToTurnOnCl:9) = nchoosek(paramArray,1); 
newStartingPoint = size(nchoosekCparamArray,1),1)+1; 
for i = 2:4 

paramisToTurnOnCnewStartingPoint:(newStartingPoint -1 + size(nchoosek(paramArray,i),1)), 
1:size(nchoosek(paramiArray,i),2))=nchoosek(paramArray,i); 
newStartingPoint = newStartingPoint + size(nchoosek(paramArray,i),1); 

end 

for j = 1:size(paramsToTurnDn,1) 

paramsToTurnOnk = paramsToTurnDn(j,:); 

y.pass the initial conditions and the parameter values to the 
yosolve_problem subroutine 
solve_problem(initcond,params); 

end 


end 

yoyosolve_problem will run and produce three plots 
yySyntax of solve_problem: solve_problem(N(0),S(0),W(0)) 
y.y.y(l) is N(t), y(2) is S(t), y(3) = W(t) 

yyoThere is a 3D plot followed by three plots in the same window, as 
yyodescribed below: 

yyoplot 1 is suppose to be the numerical solution of N 

yyoplot 2 is suppose to be the numerical solution of S 

yyoplot 3 is suppose to be the numerical solution of W 

yyoRobert Garrett Niemeyer 
yyoRichard Evan Niemeyer 

function [T,Y]=solve_problem(initCond,params) 

yyothe initial conditions are usually going to be [1,0,1] 

yyoWe set the error tolerances very low 

global a b c kmax alpha beta delta rhoO rO 
global paramsToTurnOnk 

a = params(1); 
b = params(2) ; 
c = params(3); 
kmax = params(4); 
rO = params(5); 
alpha = params(6); 
beta = params(7); 
delta = params(8); 
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rhoO = params(9); 

7,if one wants to change a, b and c at particular times, 

7othey first indicate this by changing the vector on the right 
7ohand side below to [1,2,3,0,0,0,0,0,0] 

1 

7,if one wants to change kmax rO and delta 

7othey replace the vector on the right hand side as follows to 
7othe vector [4,5,8,0,0,0,0,0,01 

7oNote, one may remove this ’if’ statement if they want to run simulate all 
7opermutations of parameter changes. 

if paramsToTurnOnk == [8,0,0,0,0,0,0,0,01 7o<- this vector here gets changed 

figure; 

options = odesetC’NonNegative’,[1 2 3],’RelTol’,le-6,’AbsTol’,le-10,’DutputFcn’,@odephas3); 
title(’test’) ; 

7oSolve the system of equations using the DDE 45 method in Matlab 

[T,Y] = ode45(@diffeqs,[1,4000],initCond,options); 

7o7oWe plot the solutions 
figure 

[theAxes,Splot,NWplot] = plotyy(T,Y(:,2),[T,Tl,[Y(:,1),Y(:,3)1 ); 
setCgca,’FontSize’,16); 

Splot.LineStyle = 

NWplot(1).LineStyle = 

NWplot(2).LineStyle = 

legend({’State resources, S ’, ’Population, N’, ’Internal conflict, W’},’FontSize’,20, . 
’Orientation’,’horizontal’,’Location’,’northoutside’); 
axes(theAxes(1)) ; 

ylabel(’State resources S’,’FontSize’,20); 
axes(theAxes(2)); 
set(gca,’FontSize’,16); 

ylabel(’Population N, Internal conflict W’,’FontSize’,20); 


end 

end 

7o7oThese are the differential equations and parameters. 
7o7oRobert Garrett Niemeyer 
7o7oRichard Evan Niemeyer 

function dy=diffeqs(t,y) 

global a b c kmax rO alpha beta delta rhoO 
global paramsToTurnOnk 



“/obelow, you can specific how you want to affect a 
“/oparticular parameter. You must also ’unaffect’ the 
“/oparEuneter below after the specification of the model 

if t>=1000 && t<2000 

if size(find(paramsToTurnOnk == 1),2) ~= 0 
a = a/3; 

end 

if size(find(paramsToTurnOnk == 2),2) ~= 0 
b = b*3; 

end 

if size(find(paramsToTurnOnk == 3),2) ~= 0 
c = c*3; 

end 

if size(find(paramsToTurnOnk == 4), 2) ~= 0 
kmax = kmax*5/3; 

end 

if size(find(paramsToTurnOnk == 5),2) ~= 0 
rO = r0*0.095/0.015; 

end 

if size(find(paramsToTurnOnk == 6),2) ~= 0 
alpha = alpha*3; 

end 

if size(find(paramsToTurnOnk == 7),2) ~= 0 
beta = beta/3; 

end 

if size(find(paramsToTurnOnk == 8),2) ~= 0 
delta = delta*9.5; 

end 

if size(find(paramsToTurnOnk == 9),2) ~= 0 
rhoO = rhoO/3; 

end 

elseif t>= 2000 && t<=4000 

if size(find(paramsToTurnOnk == 1),2) ~= 0 
a = a/9; 

end 

if size(find(paramsToTurnOnk == 2),2) ~= 0 
b = b*9; 

end 

if size(find(paramsToTurnOnk == 3),2) ~= 0 
c = c*9; 

end 

if size(find(paramsToTurnOnk == 4), 2) ~= 0 
kmax = kmax*7/3; 

end 

if size(find(paramsToTurnOnk == 5),2) ~= 0 
rO = r0*0.15/0.015; 

end 

if size(find(paramsToTurnOnk == 6),2) ~= 0 
alpha = alpha*9; 

end 

if size(find(paramsToTurnOnk == 7),2) ~= 0 
beta = beta/9; 
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end 

if size(find(paramsToTurnOnk == 8),2) ~= 0 
delta = delta*.95/.1; 

end 

if size(find(paramsToTurnOnk == 9),2) ~= 0 
rhoO = rhoO/9; 

end 

end 

dy=zeros(3,1); 


dy(l) = rO*y(l)*(l-y(l)/(kmax - c*y(3)))-delta*y(1)*y(3); 

dy(2) = rhoO*y(1)*(1-y(1)/(kmax - c*y(3)))-beta*y(l); 
dy(3) = a*y(l)~2 -b*y(3)-alpha*y(2); 

“/oBelow you will see how it is you can ’unaffect’ the parameter 
“/oby returning it to its original value. This is important. 
“/oOtherwise, the parameter value will grow exponentially larger 
“/oWith each iteration of the solver. 


if t>=1000 && t<2000 

if size(find(paramsToTurnOnk == 1),2) ~= 0 
a = a*3; 

end 

if size(find(paramsToTurnOnk == 2),2) ~= 0 
b = b/3; 

end 

if size(find(paramsToTurnOnk == 3),2) ~= 0 
c = c/3; 

end 

if size(find(paramsToTurnOnk == 4),2) ~= 0 
kmax = kmax/(5/3); 

end 

if size(find(paramsToTurnOnk == 5),2) ~= 0 
rO = r0/(0.095/0.015); 

end 

if size(find(paramsToTurnOnk == 6),2) ~= 0 
alpha = alpha/3; 

end 

if size(find(paramsToTurnOnk == 7), 2) ~= 0 
beta = beta*3; 

end 

if size(find(paramsToTurnOnk == 8),2) ~= 0 
delta = delta/9.5; 

end 

if size(find(paramsToTurnOnk == 9),2) ~= 0 
rhoO = rho0*3; 

end 

elseif t>=2000 && t<=4000 

if size(find(paramsToTurnOnk == 1),2) ~= 0 
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a = a*9; 


end 

if size(findCparamsToTurnOnk == 2),2) ~= 0 
b = b/9; 

end 

if size(findCparamsToTurnOnk == 3),2) ~= 0 
c = c/9; 

end 

if sizeCfindCparamsToTurnOnk == 4),2) ~= 0 
kmax = kmax/C7/3); 

end 

if sizeCfindCparamsToTurnOnk == 5),2) ~= 0 
rO = r0/C.15/.015); 

end 

if sizeCfindCparamsToTurnOnk == 6),2) ~= 0 
alpha = alpha/9; 

end 

if sizeCfindCparamsToTurnOnk == 7),2) ~= 0 
beta = beta*9; 

end 

if sizeCfindCparamsToTurnOnk == 8),2) ~= 0 
delta = delta/C.95/.1); 

end 

if sizeCfindCparamsToTurnOnk == 9),2) ~= 0 
rhoO = rho0*9; 

end 

end 

end 
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